# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2020/12/30 19:32
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : my_interpolation.py
@Version     : Version 1.0.0
@Description : TODO 自己实现拉格朗日插值，牛顿均差插值，埃尔米特插值，分段低次插值，三次样条插值
@Created By  : PyCharm
"""
import numpy as np
from sympy import expand, Poly, cos
from sympy.abc import x


# expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
# print(collect(expr, x))   #按x降幂排列
# print(collect(expr, x).coeff(x,2))    #和 coeff搭配使用，可以求n次项系数, 这里就是求x^2项的系数


# 自己原创
# 拉格朗日插值多项式输出显示形式与scipy.interpolate模块的同名函数有差距，有待继续优化
def lagrange(x_list: list, y_array: np.ndarray):
    """
    实现拉格朗日插值函数
    :param x_list: 插值点横坐标列表
    :param y_array: 插值点纵坐标数组
    :return: 拉格朗日插值多项式
    """
    degree = y_array.shape[0]
    l_n = 0
    for k in range(degree):
        # 注意不能同时对列表进行遍历和删除，否则索引会超出范围报错，有几种方式解决，此处使用切片
        l_n += y_array[k] * lagrange_base_function(k, x_list[::])
    return expand(l_n)


# 自己原创
def lagrange_base_function(k, x_list: list):
    """
    实现拉格朗日插值基函数
    :param k: 插值结点横坐标x_k的下标
    :param x_list: 插值点横坐标列表
    :return: 第k个插值结点的插值基函数
    """
    x_k = x_list[k]  # 取得xk处的值
    x_list.pop(k)  # 删除xk处的值
    x_list = np.array(x_list)
    #  向量运算实现插值基函数
    return np.prod((x - x_list) / (x_k - x_list))
    # 按元素实现插值基函数
    # l_k = 1
    # for i in range(x_coefficient.shape[0]):
    #     l_k *= (x - x_coefficient[i]) / (x_k - x_coefficient[i])
    # return l_k


# 自己原创
def newtown_base_function(symbol_x, x_array: np.ndarray):
    """
    实现牛顿均差插值基函数
    :param symbol_x: 符号变量x
    :param x_array: 横坐标数组
    :return: 牛顿均差插值基函数数组
    """
    nth = x_array.shape[0]
    base_function_array = np.zeros(nth, dtype=np.ndarray)
    # base_function_array = list(np.zeros(nth))
    base_function_array[0] = 1
    for i in range(1, nth):
        base_function_array[i] = np.prod(symbol_x - x_array[:i], dtype=np.ndarray)
    return base_function_array


# 自己原创
def difference_quotient_coefficient(x_array: np.ndarray, y_array: np.ndarray):
    """
    实现牛顿均差插值多项式各阶均差系数
    :param x_array: 横坐标数组
    :param y_array: 纵坐标数组
    :return: 各阶均差系数数组
    """
    nth = x_array.shape[0]
    # 初始化均差系数数组
    difference_quotient_array = np.zeros(nth)
    difference_quotient_array[0] = y_array[0]
    # 求各阶均差
    for i in range(1, nth):
        nth_difference_quotient = np.zeros(nth - i)
        # 求每阶均差
        for j in range(nth - i):
            nth_difference_quotient[j] = (y_array[j + 1] - y_array[j]) / (x_array[i + j] - x_array[j])
        difference_quotient_array[i] = nth_difference_quotient[0]
        y_array = nth_difference_quotient
    return difference_quotient_array


# 自己原创
def newtown(x_array: np.ndarray, y_array: np.ndarray, symbol_x):
    """
    实现牛顿均差插值多项式
    :param x_array: 插值结点横坐标数组
    :param y_array: 插值结点纵坐标数组
    :param symbol_x: 符号变量x
    :return: 牛顿均差插值多项式
    """
    return expand(np.sum(difference_quotient_coefficient(x_array, y_array) * newtown_base_function(symbol_x, x_array)))


# 自己原创
def qin_jiu_shao(coefficient_array: np.ndarray, value, is_ascending_order: bool = False):
    """
    多项式求值的秦九韶算法,即Horner算法
    :param coefficient_array: 多项式系数数组
    :param value: 自变量的值
    :param is_ascending_order: 多项式排列顺序
    :return: the value of polynomial
    example:
    2*x**4 - 3*x**2 + 3*x - 4=((((2x+0)x-3)x+3)x-4)
    3*x**5 - 2*x**3 + x + 7=(((((3x+0)x-2)x+0)x+1)x+7)
    """
    if is_ascending_order:
        coefficient_array = coefficient_array[::-1]
    degree = coefficient_array.shape[0]
    b = coefficient_array[0]
    for i in range(1, degree):
        b = b * value + coefficient_array[i]
    return b


if __name__ == '__main__':
    # 拉格朗日插值多项式测试成功，来源详见李庆扬数值分析第5版P48，e.x.1
    """
    x_coefficient_list = [1, -1, 2]
    x_coefficient_array = np.array(x_coefficient_list)
    y = np.array([0, -3, 4])
    lagrange_interpolation_polynomial = lagrange(x_coefficient_list, y)
    # 输出拉格朗日插值多项式
    print("拉格朗日插值多项式为:{}".format(lagrange_interpolation_polynomial))
    # 计算拉格朗日插值多项式在x=2处的值
    print("拉格朗日插值多项式在x=2处的值为:{}".format(lagrange_interpolation_polynomial.subs(x, 2)))
    """
    # 牛顿均差插值多项式测试成功，来源详见李庆扬数值分析第5版P48，e.x.1
    """
    x_coefficient_list = [1, -1, 2]
    x_coefficient_array = np.array(x_coefficient_list)
    y = np.array([0, -3, 4])
    # 牛顿均差插值多项式均差系数测试成功
    print("牛顿均差插值多项式均差系数:{}".format(difference_quotient_coefficient(x_coefficient_array, y)))
    # 牛顿均差插值多项式基函数测试成功
    print("牛顿均差插值多项式基函数:{}".format(newtown_base_function(x, x_coefficient_array)))
    newtown_polynomial = newtown(x_coefficient_array, y, x)
    # 计算牛顿插值多项式(降幂排列)在x=2处的值
    print("牛顿插值多项式(降幂排列)在x=2处的值为:{}".format(newtown_polynomial.subs(x, 2)))
    # 创建牛顿插值多项式对象
    newtown_interpolation_polynomial = Poly(newtown_polynomial, x)
    # 获得牛顿插值多项式降幂排列系数
    newtown_polynomial_coefficient = np.array(newtown_interpolation_polynomial.coeffs())  # list 类型
    # 使用秦九韶算法计算牛顿插值多项式(降幂排列)在x=2处的值
    print("秦九韶算法计算牛顿插值多项式(降幂排列)在x=2处的值为:{}".format(qin_jiu_shao(newtown_polynomial_coefficient, 2)))
    # 牛顿插值多项式次数
    print("牛顿插值多项式次数为:{}".format(newtown_interpolation_polynomial.degree()))
    # 牛顿插值多项式子项次数
    print("牛顿插值多项式子项次数为:{}".format(newtown_interpolation_polynomial.monoms()))
    """
